library(foreign)
library(readstata13)
library(effects)
library(car)
library(lattice)
library(knitr)
library(arm)
library(tidyverse)
library(dplyr)
library(plyr)
library(ggthemes)
library(tidyr)
library(Hmisc)
library(gridExtra)
library(stringr)

df <- read.dta13("coef.dta")

##########
# Figure 2
## US studies
plot.US <- ggplot(aes(x=percentile, y=coefficient), data=df) + 
  annotate(geom="text", x=50, y=50, label="Studies of the
United States",
           size=15,
           color="Black") 
plot.US + theme_few() +
  theme(panel.border = element_blank(),
        axis.title.x=element_blank(),
        axis.text.x=element_blank(),
        axis.ticks.x=element_blank(),
        axis.title.y=element_blank(),
        axis.text.y=element_blank(),
        axis.ticks.y=element_blank(),
        aspect.ratio = 1)


ggsave("US.tif", width = 180, height = 180, units = "mm", dpi=800, device="tiff")


## comparative studies 
plot.comp <- ggplot(aes(x=percentile, y=coefficient), data=df) + 
  annotate(geom="text", x=50, y=50, label="Comparative
studies",
           size=15,
           color="Black") 
plot.comp + theme_few() +
  theme(panel.border = element_blank(),
        axis.title.x=element_blank(),
        axis.text.x=element_blank(),
        axis.ticks.x=element_blank(),
        axis.title.y=element_blank(),
        axis.text.y=element_blank(),
        axis.ticks.y=element_blank(),
        aspect.ratio = 1)

ggsave("CP.tif",
       width = 180, height = 180, units = "mm", dpi=800, device="tiff")


df.Bartels2008 <- dplyr::filter(df, study=="Bartels (2008)")
plot.Bartels2008 <- ggplot(aes(x=percentile, y=coefficient), data=df.Bartels2008)+
  geom_hline(yintercept=0, color="blue", linetype="dashed", size=2) +
  geom_point(color="black", size=4, alpha = .3)+
  geom_smooth(color="black", size=3, se=FALSE)+
  scale_x_continuous(limits = c(0, 100)) +
  labs(title="Bartels (2008)", subtitle = "US, 1989-1994", x = "Income percentile", y="Coefficient")+
  theme(legend.position="FALSE") + theme_few() +
  theme(axis.text=element_text(size=16),
        axis.title=element_text(size=22),
        plot.title = element_text(size=38, face="bold", hjust=0.5),
        plot.subtitle = element_text(size=38, face="bold", hjust=0.5),
        aspect.ratio = 1)
plot.Bartels2008
ggsave("Bartels2008.tif",
       width = 180, height = 180, units = "mm", dpi=800, device="tiff")

df.Bartels2017 <- dplyr::filter(df, study=="Bartels (2017)")
plot.Bartels2017 <- ggplot(aes(x=percentile, y=coefficient), data=df.Bartels2017)+
  geom_hline(yintercept=0, color="blue", linetype="dashed", size=2) +
  geom_point(color="black", size=4, alpha = .3)+
  geom_smooth(color="black", size=3, se=FALSE)+
  scale_x_continuous(limits = c(0, 100)) +
  labs(title="Bartels (2017)", subtitle = "30 countries, 1985-2012", x = "Income percentile", y="Coefficient")+
  theme(legend.position="FALSE") + theme_few() +
  theme(axis.text=element_text(size=16),
        axis.title=element_text(size=22),
        plot.title = element_text(size=38, face="bold", hjust=0.5),
        plot.subtitle = element_text(size=38, face="bold", hjust=0.5),
        aspect.ratio = 1)
plot.Bartels2017
ggsave("Bartels2017.tif",
       width = 180, height = 180, units = "mm", dpi=800, device="tiff")

df.Bhatti2011 <- dplyr::filter(df, study=="Bhatti & Erikson (2011)")
plot.Bhatti2011 <- ggplot(aes(x=percentile, y=coefficient), data=df.Bhatti2011)+
  geom_hline(yintercept=0, color="blue", linetype="dashed", size=2) +
  geom_point(color="black", size=4, alpha = .3)+
  geom_smooth(color="black", size=3, se=FALSE, span=1.2)+
  scale_x_continuous(limits = c(0, 100)) +
  labs(title="Bhatti & Erikson (2011)", subtitle = "US, 1989-1994, 1999-2004", x = "Income percentile", y="Coefficient")+
  theme(legend.position="FALSE") + theme_few() +
  theme(axis.text=element_text(size=16),
        axis.title=element_text(size=22),
        plot.title = element_text(size=38, face="bold", hjust=0.5),
        plot.subtitle = element_text(size=38, face="bold", hjust=0.5),
        aspect.ratio = 1)
plot.Bhatti2011
ggsave("BhattiErikson2011.tif",
       width = 180, height = 180, units = "mm", dpi=800, device="tiff")


df.Brunner2013 <- dplyr::filter(df, study=="Brunner, Ross & Washington (2013)")
plot.Brunner2013 <- ggplot(aes(x=percentile, y=coefficient), data=df.Brunner2013)+
  geom_hline(yintercept=0, color="blue", linetype="dashed", size=2) +
  geom_point(color="black", size=4, alpha = .3)+
  geom_smooth(color="black", size=3, se=FALSE)+
  scale_x_continuous(limits = c(0, 100)) +
  labs(title="Brunner et al. (2013)", subtitle = "California, 1991-2008", x = "Income percentile", y="Coefficient")+
  theme(legend.position="FALSE") + theme_few() +
  theme(axis.text=element_text(size=16),
        axis.title=element_text(size=22),
        plot.title = element_text(size=38, face="bold", hjust=0.5),
        plot.subtitle = element_text(size=38, face="bold", hjust=0.5),
        aspect.ratio = 1)
plot.Brunner2013
ggsave("BrunnerRossWashington2013.tif",
       width = 180, height = 180, units = "mm", dpi=800, device="tiff")

df.Elkjaer <- dplyr::filter(df, study=="Elkjær (2020)")
plot.Elkjaer <- ggplot(aes(x=percentile, y=coefficient), data=df.Elkjaer)+
  geom_hline(yintercept=0, color="blue", linetype="dashed", size=2) +
  geom_point(color="black", size=4, alpha = .3)+
  geom_smooth(color="black", size=3, se=FALSE, span=1.8)+
  scale_x_continuous(limits = c(0, 100)) +
  labs(title="Elkjær (2020)", subtitle = "Denmark, 1985-2015", x = "Income percentile", y="Coefficient")+
  theme(legend.position="FALSE") + theme_few() +
  theme(axis.text=element_text(size=16),
        axis.title=element_text(size=22),
        plot.title = element_text(size=38, face="bold", hjust=0.5),
        plot.subtitle = element_text(size=38, face="bold", hjust=0.5),
        aspect.ratio = 1)
plot.Elkjaer
ggsave("Elkjaer2020.tif",
       width = 180, height = 180, units = "mm", dpi=800, device="tiff")


df.ElkjaerIversen <- dplyr::filter(df, study=="Elkjær & Iversen (2020)")
plot.ElkjaerIversen <- ggplot(aes(x=percentile, y=coefficient), data=df.ElkjaerIversen)+
  geom_hline(yintercept=0, color="blue", linetype="dashed", size=2) +
  geom_point(color="black", size=4, alpha = .3)+
  geom_smooth(color="black", size=3, se=FALSE)+
  scale_x_continuous(limits = c(0, 100)) +
  labs(title="Elkjær & Iversen (2020)", subtitle = "21 countries, 1985-2016", x = "Income percentile", y="Coefficient")+
  theme(legend.position="FALSE") + theme_few() +
  theme(axis.text=element_text(size=16),
        axis.title=element_text(size=22),
        plot.title = element_text(size=38, face="bold", hjust=0.5),
        plot.subtitle = element_text(size=38, face="bold", hjust=0.5),
        aspect.ratio = 1)
plot.ElkjaerIversen
ggsave("ElkjaerIversen2020.tif",
       width = 180, height = 180, units = "mm", dpi=800, device="tiff")

df.Elsasser <- dplyr::filter(df, study=="Elsässer, Hense & Schäfer (2018)")
plot.Elsasser <- ggplot(aes(x=percentile, y=coefficient), data=df.Elsasser)+
  geom_hline(yintercept=0, color="blue", linetype="dashed", size=2) +
  geom_point(color="black", size=4, alpha = .3)+
  geom_smooth(color="black", size=3, se=FALSE, span=1)+
  scale_x_continuous(limits = c(0, 100)) +
  labs(title="Elsässer et al. (2018)", subtitle = "Germany, 1998-2013", x = "Income percentile", y="Coefficient")+
  theme(legend.position="FALSE") + theme_few() +
  theme(axis.text=element_text(size=16),
        axis.title=element_text(size=22),
        plot.title = element_text(size=38, face="bold", hjust=0.5),
        plot.subtitle = element_text(size=38, face="bold", hjust=0.5),
        aspect.ratio = 1)
plot.Elsasser
ggsave("ElsasserHenseSchafer2018.tif",
       width = 180, height = 180, units = "mm", dpi=800, device="tiff")

df.Flavina <- dplyr::filter(df, study=="Flavin (2012b)")
plot.Flavina <- ggplot(aes(x=percentile, y=coefficient), data=df.Flavina)+
  geom_hline(yintercept=0, color="blue", linetype="dashed", size=2) +
  geom_point(color="black", size=4, alpha = .3)+
  geom_smooth(color="black", size=3, se=FALSE)+
  scale_x_continuous(limits = c(0, 100)) +
  labs(title="Flavin (2012a)", subtitle = "US, 1999-2004", x = "Income percentile", y="Coefficient")+
  theme(legend.position="FALSE") + theme_few() +
  theme(axis.text=element_text(size=16),
        axis.title=element_text(size=22),
        plot.title = element_text(size=38, face="bold", hjust=0.5),
        plot.subtitle = element_text(size=38, face="bold", hjust=0.5),
        aspect.ratio = 1)
plot.Flavina
ggsave("Flavin2012a.tif",
       width = 180, height = 180, units = "mm", dpi=800, device="tiff")


df.Flavinb <- dplyr::filter(df, study=="Flavin (2012a)")
plot.Flavinb <- ggplot(aes(x=percentile, y=coefficient), data=df.Flavinb)+
  geom_hline(yintercept=0, color="blue", linetype="dashed", size=2) +
  geom_point(color="black", size=4, alpha = .3)+
  geom_smooth(color="black", size=3, se=FALSE)+
  scale_x_continuous(limits = c(0, 100)) +
  labs(title="Flavin (2012b)", subtitle = "US, 1988-1992, 2000-2004", x = "Income percentile", y="Coefficient")+
  theme(legend.position="FALSE") + theme_few() +
  theme(axis.text=element_text(size=16),
        axis.title=element_text(size=22),
        plot.title = element_text(size=38, face="bold", hjust=0.5),
        plot.subtitle = element_text(size=38, face="bold", hjust=0.5),
        aspect.ratio = 1)
plot.Flavinb
ggsave("Flavin2012b.tif",
       width = 180, height = 180, units = "mm", dpi=800, device="tiff")

df.Gilens2005 <- dplyr::filter(df, study=="Gilens (2005)")
plot.Gilens2005 <- ggplot(aes(x=percentile, y=coefficient), data=df.Gilens2005)+
  geom_hline(yintercept=0, color="blue", linetype="dashed", size=2) +
  geom_point(color="black", size=4, alpha = .3)+
  geom_smooth(color="black", size=3, se=FALSE, span=1.2)+
  scale_x_continuous(limits = c(0, 100)) +
  labs(title="Gilens (2005)", subtitle = "US, 1981-2002", x = "Income percentile", y="Coefficient")+
  theme(legend.position="FALSE") + theme_few() +
  theme(axis.text=element_text(size=16),
        axis.title=element_text(size=22),
        plot.title = element_text(size=38, face="bold", hjust=0.5),
        plot.subtitle = element_text(size=38, face="bold", hjust=0.5),
        aspect.ratio = 1)
plot.Gilens2005
ggsave("Gilens2005.tif",
       width = 180, height = 180, units = "mm", dpi=800, device="tiff")

df.Gilens2012 <- dplyr::filter(df, study=="Gilens (2012)")
plot.Gilens2012 <- ggplot(aes(x=percentile, y=coefficient), data=df.Gilens2012)+
  geom_hline(yintercept=0, color="blue", linetype="dashed", size=2) +
  geom_point(color="black", size=4, alpha = .3)+
  geom_smooth(color="black", size=3, se=FALSE)+
  scale_x_continuous(limits = c(0, 100)) +
  labs(title="Gilens (2012)", subtitle = "US, 1964-1968, 1981-2002", x = "Income percentile", y="Coefficient")+
  theme(legend.position="FALSE") + theme_few() +
  theme(axis.text=element_text(size=16),
        axis.title=element_text(size=22),
        plot.title = element_text(size=38,face="bold", hjust=0.5),
        plot.subtitle = element_text(size=38, face="bold",hjust=0.5),
        aspect.ratio = 1)
plot.Gilens2012
ggsave("Gilens2012.tif",
       width = 180, height = 180, units = "mm", dpi=800, device="tiff")

df.Hayes <- dplyr::filter(df, study=="Hayes (2012)")
plot.Hayes <- ggplot(aes(x=percentile, y=coefficient), data=df.Hayes)+
  geom_hline(yintercept=0, color="blue", linetype="dashed", size=2) +
  geom_point(color="black", size=4, alpha = .3)+
  geom_smooth(color="black", size=3, se=FALSE)+
  scale_x_continuous(limits = c(0, 100)) +
  labs(title="Hayes (2012)", subtitle = "US, 2001-2011", x = "Income percentile", y="Coefficient")+
  theme(legend.position="FALSE") + theme_few() +
  theme(axis.text=element_text(size=16),
        axis.title=element_text(size=22),
        plot.title = element_text(size=38, face="bold", hjust=0.5),
        plot.subtitle = element_text(size=38, face="bold", hjust=0.5),
        aspect.ratio = 1)
plot.Hayes
ggsave("Hayes2012.tif",
       width = 180, height = 180, units = "mm", dpi=800, device="tiff")

df.Peters <- dplyr::filter(df, study=="Peters & Ensink (2015)")
plot.Peters <- ggplot(aes(x=percentile, y=coefficient), data=df.Peters)+
  geom_hline(yintercept=0, color="blue", linetype="dashed", size=2) +
  geom_point(color="black", size=4, alpha = .3)+
  geom_smooth(color="black", size=3, se=FALSE)+
  scale_x_continuous(limits = c(0, 100)) +
  labs(title="Peters & Ensink (2015)", subtitle = "25 countries, 2002-2010", x = "Income percentile", y="Coefficient")+
  theme(legend.position="FALSE") + theme_few() +
  theme(axis.text=element_text(size=16),
        axis.title=element_text(size=22),
        plot.title = element_text(size=38, face="bold", hjust=0.5),
        plot.subtitle = element_text(size=38, face="bold", hjust=0.5),
        aspect.ratio = 1)
plot.Peters
ggsave("PetersEnsink2015.tif",
       width = 180, height = 180, units = "mm", dpi=800, device="tiff")


df.Rigby2011 <- dplyr::filter(df, study=="Rigby & Wright (2011)")
plot.Rigby2011 <- ggplot(aes(x=percentile, y=coefficient), data=df.Rigby2011)+
  geom_hline(yintercept=0, color="blue", linetype="dashed", size=2) +
  geom_point(color="black", size=4, alpha = .3)+
  geom_smooth(color="black", size=3, se=FALSE, span=1.1)+
  scale_x_continuous(limits = c(0, 100)) +
  labs(title="Rigby & Wright (2011)", subtitle = "US, 2000-2004", x = "Income percentile", y="Coefficient")+
  theme(legend.position="FALSE") + theme_few() +
  theme(axis.text=element_text(size=16),
        axis.title=element_text(size=22),
        plot.title = element_text(size=38, face="bold", hjust=0.5),
        plot.subtitle = element_text(size=38, face="bold", hjust=0.5),
        aspect.ratio = 1)
plot.Rigby2011
ggsave("RigbyWright2011.tif",
       width = 180, height = 180, units = "mm", dpi=800, device="tiff")


df.Rigby2013 <- dplyr::filter(df, study=="Rigby & Wright (2013)")
plot.Rigby2013 <- ggplot(aes(x=percentile, y=coefficient), data=df.Rigby2013)+
  geom_hline(yintercept=0, color="blue", linetype="dashed", size=2) +
  geom_point(color="black", size=4, alpha = .3)+
  geom_smooth(color="black", size=3, se=FALSE)+
  scale_x_continuous(limits = c(0, 100)) +
  labs(title="Rigby & Wright (2013)", subtitle = "US, 2000", x = "Income percentile", y="Coefficient")+
  theme(legend.position="FALSE") + theme_few() +
  theme(axis.text=element_text(size=16),
        axis.title=element_text(size=22),
        plot.title = element_text(size=38, face="bold", hjust=0.5),
        plot.subtitle = element_text(size=38, face="bold", hjust=0.5),
        aspect.ratio = 1)
plot.Rigby2013
ggsave("RigbyWright2013.tif",
       width = 180, height = 180, units = "mm", dpi=800, device="tiff")


df.Soroka <- dplyr::filter(df, study=="Soroka & Wlezien (2010)")
plot.Soroka <- ggplot(aes(x=percentile, y=coefficient), data=df.Soroka)+
  geom_hline(yintercept=0, color="blue", linetype="dashed", size=2) +
  geom_point(color="black", size=4, alpha = .3)+
  geom_smooth(color="black", size=3, se=FALSE)+
  scale_x_continuous(limits = c(0, 100)) +
  labs(title="Soroka & Wlezien (2010)", subtitle = "US & Canada, 1973-2005", x = "Income percentile", y="Coefficient")+
  theme(legend.position="FALSE") + theme_few() +
  theme(axis.text=element_text(size=16),
        axis.title=element_text(size=22),
        plot.title = element_text(size=38, face="bold", hjust=0.5),
        plot.subtitle = element_text(size=38, face="bold", hjust=0.5),
        aspect.ratio = 1)
plot.Soroka
ggsave("SorokaWlezien2010.tif",
       width = 180, height = 180, units = "mm", dpi=800, device="tiff")

df.Wlezien <- dplyr::filter(df, study=="Wlezien & Soroka (2011)")
plot.Wlezien <- ggplot(aes(x=percentile, y=coefficient), data=df.Wlezien)+
  geom_hline(yintercept=0, color="blue", linetype="dashed", size=2) +
  geom_point(color="black", size=4, alpha = .3)+
  geom_smooth(color="black", size=3, se=FALSE)+
  scale_x_continuous(limits = c(0, 100)) +
  labs(title="Wlezien & Soroka (2011)", subtitle = "US, 1973-2008", x = "Income percentile", y="Coefficient")+
  theme(legend.position="FALSE") + theme_few() +
  theme(axis.text=element_text(size=16),
        axis.title=element_text(size=22),
        plot.title = element_text(size=38, face="bold", hjust=0.5),
        plot.subtitle = element_text(size=38, face="bold", hjust=0.5),
        aspect.ratio = 1)
plot.Wlezien
ggsave("WlezienSoroka2011.tif",
       width = 180, height = 180, units = "mm", dpi=800, device="tiff")

df.Tausanovitch <- dplyr::filter(df, study=="Tausanovitch (2016)")
plot.Tausanovitch <- ggplot(aes(x=percentile, y=coefficient), data=df.Tausanovitch)+
  geom_hline(yintercept=0, color="blue", linetype="dashed", size=2) +
  geom_point(color="black", size=4, alpha = .3)+
  geom_smooth(color="black", size=3, se=FALSE)+
  scale_x_continuous(limits = c(0, 100)) +
  labs(title="Tausanovitch (2016)", subtitle = "US, 2000-2012", x = "Income percentile", y="Coefficient")+
  theme(legend.position="FALSE") + theme_few() +
  theme(axis.text=element_text(size=16),
        axis.title=element_text(size=22),
        plot.title = element_text(size=38, face="bold", hjust=0.5),
        plot.subtitle = element_text(size=38, face="bold", hjust=0.5),
        aspect.ratio = 1)
plot.Tausanovitch
ggsave("Tausanovitch2016.tif",
       width = 180, height = 180, units = "mm", dpi=800, device="tiff")

df.Ura <- dplyr::filter(df, study=="Ura & Ellis (2008)")
plot.Ura <- ggplot(aes(x=percentile, y=coefficient), data=df.Ura)+
  geom_hline(yintercept=0, color="blue", linetype="dashed", size=2) +
  geom_point(color="black", size=4, alpha = .3)+
  geom_smooth(color="black", size=3, se=FALSE)+
  scale_x_continuous(limits = c(0, 100)) +
  labs(title="Ura & Ellis (2008)", subtitle = "US, 1974-1996", x = "Income percentile", y="Coefficient")+
  theme(legend.position="FALSE") + theme_few() +
  theme(axis.text=element_text(size=16),
        axis.title=element_text(size=22),
        plot.title = element_text(size=38, face="bold", hjust=0.5),
        plot.subtitle = element_text(size=38, face="bold", hjust=0.5),
        aspect.ratio = 1)
plot.Ura
ggsave("UraEllis2008.tif",
       width = 180, height = 180, units = "mm", dpi=800, device="tiff")


df.GilensPage2014 <- dplyr::filter(df, study=="Gilens & Page (2014)")
plot.GilensPage2014 <- ggplot(aes(x=percentile, y=coefficient), data=df.GilensPage2014)+
  geom_hline(yintercept=0, color="blue", linetype="dashed", size=2) +
  geom_point(color="black", size=4, alpha = .3)+
  geom_smooth(color="black", size=3, se=FALSE)+
  scale_x_continuous(limits = c(0, 100)) +
  labs(title="Gilens & Page (2014)", subtitle = "US, 1981-2002", x = "Income percentile", y="Coefficient")+
  theme(legend.position="FALSE") + theme_few() +
  theme(axis.text=element_text(size=16),
        axis.title=element_text(size=22),
        plot.title = element_text(size=38, face="bold", hjust=0.5),
        plot.subtitle = element_text(size=38, face="bold", hjust=0.5),
        aspect.ratio = 1)
plot.GilensPage2014
ggsave("GilensPage2014.tif",
       width = 180, height = 180, units = "mm", dpi=800, device="tiff")

df.Rhodes_Schaffner <- dplyr::filter(df, study=="Rhodes & Schaffner (2017)")
plot.Rhodes_Schaffner <- ggplot(aes(x=percentile, y=coefficient), data=df.Rhodes_Schaffner)+
  geom_hline(yintercept=0, color="blue", linetype="dashed", size=2) +
  geom_point(color="black", size=4, alpha = .3)+
  geom_smooth(color="black", size=3, se=FALSE)+
  scale_x_continuous(limits = c(0, 100)) +
  labs(title="Rhodes & Schaffner (2017)", subtitle = "US, 2012", x = "Income percentile", y="Coefficient")+
  theme(legend.position="FALSE") + theme_few() +
  theme(axis.text=element_text(size=16),
        axis.title=element_text(size=22),
        plot.title = element_text(size=38, face="bold", hjust=0.6),
        plot.subtitle = element_text(size=38, face="bold", hjust=0.5),
        aspect.ratio = 1)
plot.Rhodes_Schaffner
ggsave("RhodesSchaffner.tif",
       width = 180, height = 180, units = "mm", dpi=800, device="tiff")


df.Schakel <- dplyr::filter(df, study=="Schakel (2019)")
plot.Schakel <- ggplot(aes(x=percentile, y=coefficient), data=df.Schakel)+
  geom_hline(yintercept=0, color="blue", linetype="dashed", size=2) +
  geom_point(color="black", size=4, alpha = .3)+
  geom_smooth(color="black", size=3, se=FALSE)+
  scale_x_continuous(limits = c(0, 100)) +
  labs(title="Schakel (2019)", subtitle = "Netherlands, 1979-2012", x = "Income percentile", y="Coefficient")+
  theme(legend.position="FALSE") + theme_few() +
  theme(axis.text=element_text(size=16),
        axis.title=element_text(size=22),
        plot.title = element_text(size=38, face="bold", hjust=0.5),
        plot.subtitle = element_text(size=38, face="bold", hjust=0.5),
        aspect.ratio = 1)
plot.Schakel
ggsave("Schakel2019.tif",
       width = 180, height = 180, units = "mm", dpi=800, device="tiff")


df.SchakelBurgoonHakhverdian <- dplyr::filter(df, study=="Schakel, Burgoon & Hakhverdian (2020)")
plot.SchakelBurgoonHakhverdian <- ggplot(aes(x=percentile, y=coefficient), data=df.SchakelBurgoonHakhverdian)+
  geom_hline(yintercept=0, color="blue", linetype="dashed", size=2) +
  geom_point(color="black", size=4, alpha = .3)+
  geom_smooth(color="black", size=3, se=FALSE)+
  scale_x_continuous(limits = c(0, 100)) +
  labs(title="Schakel et al. (2020)", subtitle = "20 countries, 1985-2006", x = "Income percentile", y="Coefficient")+
  theme(legend.position="FALSE") + theme_few() +
  theme(axis.text=element_text(size=16),
        axis.title=element_text(size=22),
        plot.title = element_text(size=38, face="bold", hjust=0.5),
        plot.subtitle = element_text(size=38, face="bold", hjust=0.5),
        aspect.ratio = 1)
plot.SchakelBurgoonHakhverdian
ggsave("SchakelBurgoonHakhverdian2020.tif",
       width = 180, height = 180, units = "mm", dpi=800, device="tiff")

df.Lax <- dplyr::filter(df, study=="Lax, Phillips & Zelizer (2019)")
plot.Lax <- ggplot(aes(x=percentile, y=coefficient), data=df.Lax)+
  geom_hline(yintercept=0, color="blue", linetype="dashed", size=2) +
  geom_point(color="black", size=4, alpha = .3)+
  geom_smooth(color="black", size=3, se=FALSE)+
  scale_x_continuous(limits = c(0, 100)) +
  labs(title="Lax et al. (2019)", subtitle = "US, 2001-2015", x = "Income percentile", y="Coefficient")+
  theme(legend.position="FALSE") + theme_few() +
  theme(axis.text=element_text(size=16),
        axis.title=element_text(size=22),
        plot.title = element_text(size=38, face="bold", hjust=0.5),
        plot.subtitle = element_text(size=38, face="bold", hjust=0.5),
        aspect.ratio = 1)
plot.Lax
ggsave("LaxPhillipsZelizer2019.tif",
       width = 180, height = 180, units = "mm", dpi=800, device="tiff")


df.Stadelmann <- dplyr::filter(df, study=="Stadelmann, Portmann & Eichenberger (2015)")
plot.Stadelmann <- ggplot(aes(x=percentile, y=coefficient), data=df.Stadelmann)+
  geom_hline(yintercept=0, color="blue", linetype="dashed", size=2) +
  geom_point(color="black", size=4, alpha = .3)+
  geom_smooth(color="black", size=3, se=FALSE)+
  scale_x_continuous(limits = c(0, 100)) +
  labs(title="Stadelmann et al. (2015)", subtitle = "Switzerland, 1996-2012", x = "Income percentile", y="Coefficient")+
  theme(legend.position="FALSE") + theme_few() +
  theme(axis.text=element_text(size=16),
        axis.title=element_text(size=22),
        plot.title = element_text(size=38, face="bold", hjust=0.5),
        plot.subtitle = element_text(size=38, face="bold", hjust=0.5),
        aspect.ratio = 1)
plot.Stadelmann
ggsave("StadelmannPortmannEichenberger2015.tif",
       width = 180, height = 180, units = "mm", dpi=800, device="tiff")

df.Wright <- dplyr::filter(df, study=="Wright & Rigby (2020)")
plot.Wright <- ggplot(aes(x=percentile, y=coefficient), data=df.Wright)+
  geom_hline(yintercept=0, color="blue", linetype="dashed", size=2) +
  geom_point(color="black", size=4, alpha = .3)+
  geom_smooth(color="black", size=3, se=FALSE)+
  scale_x_continuous(limits = c(0, 100)) +
  labs(title="Wright & Rigby  (2020)", subtitle = "US, 2000-2004", x = "Income percentile", y="Coefficient")+
  theme(legend.position="FALSE") + theme_few() +
  theme(axis.text=element_text(size=16),
        axis.title=element_text(size=22),
        plot.title = element_text(size=38, face="bold", hjust=0.5),
        plot.subtitle = element_text(size=38, face="bold", hjust=0.5),
        aspect.ratio = 1)
plot.Wright
ggsave("WrightRigby.tif",
       width = 180, height = 180, units = "mm", dpi=800, device="tiff")

####################################################################################################
# Figure 3
df.wide <- read.dta13("coef_wide.dta")
df.wide2 <- dplyr::filter(df.wide, t_L<10 & t_M<10 & t_H<10)

p.dens <- ggplot(data=df.wide2) +
  geom_histogram(aes(x =t_L), bins=50, linetype=0, colour="red", fill="red", alpha=.3) +
  geom_histogram(aes(x =t_M), bins=50, linetype=0, colour="darkgreen", fill="darkgreen", alpha=.3) +
  geom_histogram(aes(x =t_H), bins=50, linetype=0, colour="blue", fill="blue", alpha=.3) +
  scale_x_continuous(breaks = c(-5,-1.96,0,1.96,5,10,15), labels = c("-5","-1.96","0","1.96","5","10","15")) +
  labs(x = expression(paste("Test statistic (", beta/se,")")), y="")+
  theme_bw() +
  theme(axis.text=element_text(size=16),
        axis.title=element_text(size=18))

## text for figure 3
p.dens +
  annotate(geom="text", x=-3, y=14.25, label=expression("low-income"),
           size=7,
           color="red", 
           alpha=.9) +
  annotate(geom="text", x=-3, y=13.5, label=expression("test statistics"),
           size=7,
           color="red", 
           alpha=.9) +
  geom_curve(aes(x=-2.65,y=13,xend=-1.5,yend=11.5),
             curvature = 0.1, linetype="solid", size=.1, color="red",
             arrow = arrow(type="open",length = unit(0.25,"cm"))) +
  annotate(geom="text", x=1.4, y=17.25, label=expression("middle-income"),
           size=7,
           color="darkgreen", 
           alpha=.9) +
  annotate(geom="text", x=1.4, y=16.5, label=expression("test statistics"),
           size=7,
           color="darkgreen", 
           alpha=.9) +
  geom_curve(aes(x=1.3,y=15.6,xend=.85,yend=13.25),
             curvature = 0.1, linetype="solid", size=.1, color="darkgreen",
             arrow = arrow(type="open",length = unit(0.25,"cm"))) +
  annotate(geom="text", x=5.5, y=13, label=expression("high-income"),
           size=7,
           color="blue", 
           alpha=.9) +
  annotate(geom="text", x=5.5, y=12.25, label=expression("test statistics"),
           size=7,
           color="blue", 
           alpha=.9) +
  geom_curve(aes(x=4.75,y=11.75,xend=3.75,yend=10.5),
             curvature = -0.1, linetype="solid", size=.1, color="blue",
             arrow = arrow(type="open",length = unit(0.25,"cm"))) 


ggsave("t_statistics.tif",
       width = 250, height = 200, units = "mm", dpi=800, device="tiff")


####################################################################################################
# text for figure 4
df.text.A <- data.frame(
  label = "A",
  x     = .75,
  y     = .32
)
df.text.B <- data.frame(
  label = "B",
  x     = .75,
  y     = .32
)

# Figure 4a
p.HL <- ggplot(data = df.wide) + 
  geom_bar(mapping = aes(x = HL_7_sig, y = ..prop.., HL_7_sig = 1), 
           stat = "count", width = .9, color="black", 
           fill=c("#5F3B3E", "#BF777D", "#ef959d", "#fcddbc","#cfdbbd" , "#b8d8ba", "#9cac9d", "#a59cac"),
           size=1, 
           alpha=1)+
  geom_vline(xintercept=7.5, color="red", linetype="dashed", size=1) +
  ylab("Proportion") +
  xlab("") +
  scale_y_continuous(limits = c(0, .35),
                     expand = c(0.01, 0)) +
  scale_x_continuous(guide = guide_axis(n.dodge = 1), 
                     labels = function(x) stringr::str_wrap(c("extreme pro-poor bias", 
                                                              "severe pro-poor bias", 
                                                              "moderate pro-poor bias",
                                                              "relatively equal representation", 
                                                              "moderate pro-rich bias", 
                                                              "severe pro-rich bias",
                                                              "extreme pro-rich bias",
                                                              "ambiguous results"), width = 10),
                     breaks=c(1,2,3,4,5,6,7,8))+
  geom_text(
    data    = df.text.A,
    mapping = aes(x = x, y = y, label = label), size=15, fontface="bold"
  )+
  theme_bw() +
  theme(axis.text=element_text(size=20),
        axis.title=element_text(size=24),
        axis.ticks = element_blank())

p.HL

ggsave("HL.tif",
       width = 400, height = 125, units = "mm", dpi=800, device="tiff")


# Figure 4b
p.HM <- ggplot(data = df.wide) + 
  geom_bar(mapping = aes(x = HM_7_sig, y = ..prop.., HM_7_sig = 1), 
           stat = "count", width = .9, color="black", 
           fill=c("#5F3B3E", "#BF777D", "#ef959d", "#fcddbc","#cfdbbd" , "#b8d8ba", "#9cac9d", "#a59cac"),
           size=1, 
           alpha=1)+
  geom_vline(xintercept=7.5, color="red", linetype="dashed", size=1) +
  ylab("Proportion") +
  xlab("") +
  scale_y_continuous(limits = c(0, .35),
                     expand = c(0.01, 0)) +
  scale_x_continuous(guide = guide_axis(n.dodge = 1), 
                     labels = function(x) stringr::str_wrap(c("extreme pro-middle-class bias", 
                                                              "severe pro-middle-class bias", 
                                                              "moderate pro-middle-class bias",
                                                              "relatively equal representation", 
                                                              "moderate pro-rich bias", 
                                                              "severe pro-rich bias",
                                                              "extreme pro-rich bias",
                                                              "ambiguous results"), width = 10),
                     breaks=c(1,2,3,4,5,6,7,8))+
  geom_text(
    data    = df.text.B,
    mapping = aes(x = x, y = y, label = label), size=15, fontface="bold"
  )+
  theme_bw() +
  theme(axis.text=element_text(size=20),
        axis.title=element_text(size=24),
        axis.ticks = element_blank())

p.HM
ggsave("HM.tif",
       width = 400, height = 125, units = "mm", dpi=800, device="tiff")

#### Figures 5 and 6: Predicted probabilities
fill_color <- c("#9cac9d", "#b8d8ba","#d9dbbc","#fcddbc", "#ef959d", "#A9888B","#69585f")


### H-L



#### Figures 5-9: Predicted probabilities
fill_color <- c("#9cac9d", "#b8d8ba","#cfdbbd","#fcddbc", "#ef959d", "#BF777D","#5F3B3E")
### H-M/L
HL1 <- "extreme pro-poor/ \n middle-class bias" 
HL2 <- "severe pro-poor/ \n middle-class bias"
HL3 <- "moderate pro-poor/ \n middle-class bias"
HL4 <- "relatively equal \n representation"
HL5 <- "moderate pro-rich bias"
HL6 <- "severe pro-rich bias"
HL7 <- "extreme pro-rich bias"

# Region
df.all <- read.dta13("pr_comb.dta")
df.region <- dplyr::filter(df.all, byvar=="Region")

p <- ggplot(df.region, aes(x = re_lab, y = pr, fill = as.factor(-HL_7))) + 
  geom_bar(position="fill", stat = "identity")+ 
  facet_wrap(.~class, nrow=1)+
  scale_fill_manual(values= fill_color, guide = guide_legend(reverse=FALSE, nrow = 7), name = "", labels = c(HL7, HL6, HL5, HL4, HL3, HL2, HL1)) +
  scale_y_continuous(breaks = c(0,.2,.4,.6,.8,1), limits = c(0, 1),
                     expand = c(0, 0)) +
  labs(title="", x = "Region", y="Predicted probability", fill="")+
  theme_bw() +
  theme(axis.text=element_text(size=15),
        axis.title=element_text(size=16),
        strip.background = element_rect(colour="white", fill="white"),
        legend.position = "right",
        legend.text = element_text(size=15),
        strip.text.x = element_text(size = 16, face="bold"),
        legend.key.size = unit(3.1, 'lines'))
p 

ggsave("region.tif",
       width = 320, height = 180, units = "mm", dpi=800, device="tiff")

# N of groups
df.groups <- dplyr::filter(df.all, byvar=="N of groups included in a model")
df.groups.sorted <- mutate(df.groups, re_lab = reorder(re_lab, lab))


p <- ggplot(df.groups.sorted, aes(x = re_lab, y = pr, fill = as.factor(-HL_7))) + 
  geom_bar(position="fill", stat = "identity")+ 
  facet_wrap(.~class, nrow=1)+
  scale_fill_manual(values= fill_color, guide = guide_legend(reverse=FALSE, nrow = 7), name = "", labels = c(HL7, HL6, HL5, HL4, HL3, HL2, HL1)) +
  scale_y_continuous(breaks = c(0,.2,.4,.6,.8,1), limits = c(0, 1),
                     expand = c(0, 0)) +
  labs(title="", x = "Number of groups included in a model", y="Predicted probability", fill="")+
  theme_bw() +
  theme(axis.text=element_text(size=15),
        axis.title=element_text(size=16),
        strip.background = element_rect(colour="white", fill="white"),
        legend.position = "right",
        legend.text = element_text(size=15),
        strip.text.x = element_text(size = 16, face="bold"),
        legend.key.size = unit(3.1, 'lines'))
p 

ggsave("N_groups.tif",
       width = 320, height = 180, units = "mm", dpi=800, device="tiff")

# Partisanship
df.party <- dplyr::filter(df.all, byvar=="Partisanship")
df.party.sorted <- mutate(df.party, re_lab = reorder(re_lab, lab))


p <- ggplot(df.party.sorted, aes(x = re_lab, y = pr, fill = as.factor(-HL_7))) + 
  geom_bar(position="fill", stat = "identity")+ 
  facet_wrap(.~class, nrow=1)+
  scale_fill_manual(values= fill_color, guide = guide_legend(reverse=FALSE, nrow = 7), name = "", labels = c(HL7, HL6, HL5, HL4, HL3, HL2, HL1)) +
  scale_y_continuous(breaks = c(0,.2,.4,.6,.8,1), limits = c(0, 1),
                     expand = c(0, 0)) +
  labs(title="", x = "Partisanship", y="Predicted probability", fill="")+
  theme_bw() +
  theme(axis.text=element_text(size=15),
        axis.title=element_text(size=16),
        strip.background = element_rect(colour="white", fill="white"),
        legend.position = "right",
        legend.text = element_text(size=15),
        strip.text.x = element_text(size = 16, face="bold"),
        legend.key.size = unit(3.1, 'lines'))
p 

ggsave("party.tif",
       width = 320, height = 180, units = "mm", dpi=800, device="tiff")

# Policy domain
df.issue <- dplyr::filter(df.all, byvar=="Issue")
df.issue.sorted <- mutate(df.issue, re_lab = reorder(re_lab, lab))

p <- ggplot(df.issue.sorted, aes(x = re_lab, y = pr, fill = as.factor(-HL_7))) + 
  geom_bar(position="fill", stat = "identity")+ 
  facet_wrap(.~class, nrow=1)+
  scale_fill_manual(values= fill_color, guide = guide_legend(reverse=FALSE, nrow = 7), name = "", labels = c(HL7, HL6, HL5, HL4, HL3, HL2, HL1)) +
  scale_y_continuous(breaks = c(0,.2,.4,.6,.8,1), limits = c(0, 1),
                     expand = c(0, 0)) +
  labs(title="", x = "Policy domain", y="Predicted probability", fill="")+
  theme_bw() +
  theme(axis.text=element_text(size=15),
        axis.title=element_text(size=16),
        strip.background = element_rect(colour="white", fill="white"),
        legend.position = "right",
        legend.text = element_text(size=15),
        strip.text.x = element_text(size = 16, face="bold"),
        legend.key.size = unit(3.1, 'lines'))
p 

ggsave("issue.tif",
       width = 320, height = 180, units = "mm", dpi=800, device="tiff")


# General Political Ideology
df.GPI <- dplyr::filter(df.all, byvar=="General political ideology")
df.GPI.sorted <- mutate(df.GPI, re_lab = reorder(re_lab, lab))

p <- ggplot(df.GPI.sorted, aes(x = re_lab, y = pr, fill = as.factor(-HL_7))) + 
  geom_bar(position="fill", stat = "identity")+ 
  facet_wrap(.~class, nrow=1)+
  scale_fill_manual(values= fill_color, guide = guide_legend(reverse=FALSE, nrow = 7), name = "", labels = c(HL7, HL6, HL5, HL4, HL3, HL2, HL1)) +
  scale_y_continuous(breaks = c(0,.2,.4,.6,.8,1), limits = c(0, 1),
                     expand = c(0, 0)) +
  labs(title="", x = "General political ideology", y="Predicted probability", fill="")+
  theme_bw() +
  theme(axis.text=element_text(size=15),
        axis.title=element_text(size=16),
        strip.background = element_rect(colour="white", fill="white"),
        legend.position = "right",
        legend.text = element_text(size=15),
        strip.text.x = element_text(size = 16, face="bold"),
        legend.key.size = unit(3.1, 'lines'))
p 
ggsave("global.tif",
       width = 320, height = 180, units = "mm", dpi=800, device="tiff")
